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>~>' Abstract 

^. 

I '^ ' In this paper we investigate a subgrid model based on an anisotropic version 

,jj. _ of the NS-a model using a lid-driven cavity flow at a Reynolds number of 10,000. 

^ I Previously the NS-a model has only been used numerically in the isotropic form. 

'NT . 

lO . The subgrid model is developed from the Eulerian-averaged anisotropic equations 

f^ ■ [Holm, Physica D, v. 133, pp 215-269, 1999]. It was found that when a^ was based 

on the mesh numerical oscillations developed which manifested themselves in the 
[~^ I appearance of streamwise vortices and a 'mixing out' of the velocity profile. This 

is analogous to the Craik-Leibovich mechanism, with the difference being that the 



% 



oscillations here are not physical but numerical. The problem could be traced back 
to the discontinuity in a^ encountered when a^ = on the endwalls. An alternative 
definition of a^ based on velocity gradients, rather than mesh spacing, is proposed 
and tested. Using this definition the results with the model shown a significant 
improvement. The splitting of the downstream wall jet, rms and shear stress profiles 
are correctly captured a coarse mesh. The model is shown to predict both positive 
and negative energy transfer in the jet impingement region, in qualitative agreement 
with DNS results. 



1 Introduction 

An accurate description of turbulent flows is of paramount importance both in terms of 
engineering applications, and in understanding physical phenomena in the natural world. 
Increasingly, numerical computations are playing a prominent role in turbulence research. 
However, for many practical problems, the full range of scales active in a turbulent flow 
cannot be resolved on a finite computational domain. This means models must be intro- 
duced to parameterize the effects of the unresolved motions on the resolved ones. This is 
usually done by first applying averaging procedures directly to the Navier-Stokes equa- 
tions. The most common methods are either to introduce a statistical average, which 
leads to the Reynolds-Averaged Navier-Stokes (RANS) equations, or to use a spatial fil- 
ter, which leads to the Large-Eddy Simulation (LES) equations. Both methods lead to 
the appearance of an unclosed term, for which a wide variety of models have been pre- 
sented for both RANS and LES methodologies. ^^'^'^ One consistent trend has been the 
use of models which are either strictly dissipative, or contain a dissipative component. 
This concept is well-founded, since the role of the small scales, which are being modeled, 
is to remove the energy generated through non-linear interactions of the large, resolved 
scales. However, there are some fiows where the non-linear interactions are weak, and the 
dissipation provided by such models may be excessive. An example of this is the early 
stages of transition in a boundary layer fiow, where dissipative models may delay, or even 
prevent, the onset of transition. ^^ 

The NS-a model is considered a non-dissiptive model in that it uses a modified nonlinear- 
ity to alter the energy transfer among scales instead of changing the dissipative process. ^^ 
However, it has different origins than other nonlinear models. Instead of starting with the 
Navier-Stokes equations, the governing equations can be derived by applying Hamilton's 
principle to an averaged Lagrangian,^^ and the resulting (inviscid) equations conserve 
circulation when the model parameter o? is constant. The derivation using Hamilton's 
principle naturally leads to a set of equations which contain two velocity fields ui and -Uj, 
where ui is smoother than ui via an inversion of the Helmholtz operator. In the isotropic 
case the the parameter that arises in the averaging procedure is a scalar, cP' . When c^ is 
constant the governing equations can be written^ 

diUi = 0, (1) 



dtUi + UjdjUi + UkdiUk = -dip°' + udkkUi, (2) 

with, 

Ui= {l- a^dkk) Ui, (3) 

1 
p"" =p- -UiUi. (4) 

The third term on the LHS of (2) is unique to the NS-a model and will be referred to in 
the following as the tilting term because it arises due to the velocity difference between two 
ends of a Lagrangian trajectory which is being carried by a smoothed flow. We also use 
the standard Laplacian operator acting on the momentum velocity Ui in the dissipation 
term in the interest of maintaining a model similar to that used in other studies.^' ^^'^^'^^ 
The above set of equations are also known as the viscous Camassa-Holm equations. In- 
terest in using these equations as a model for turbulence can be traced back to the work 
by Chen^'^ where analytical results shown to yield velocity and shear stress profiles in 
good agreement with experimental results for pipe and channel flows. 

The early literature on the NS-a equations hypothesized that the equations would have 
an energy spectrum with a steeper slope in the inertial subrange for length scales smaller 
than a,*^'^'^'^ making it a good candidate for an LES model. The slope based on the 
conserved energy E^ = jyUiUidV was expected to be k~^, which corresponds to k~^ for 
the translational energy E^ = jy u-iUi dV . The physical mechanism behind the steeper 
slope was explained as the suppression of nonlinear interactions between scales which are 
smaller than a.^ More recently it has been shown in an enlightening study^^ that the en- 
ergy spectrum is in fact not steeper than that of the Navier-Stokes equations. The reason 
for this is that the NS-a fluid is comprised of both regions undergoing the Navier-Stokes 
dynamics of vorticity transport and stretching, and regions described as 'rigid rotators'^'* 
where stretching is inhibited. These rigid rotators have no internal degrees of freedom 
but do have kinetic energy. Their scaling leads to an energy spectrum which has a slope 
of /c^, which means the k~^ spectrum is subdominant.^^ Thus while the hypothesis of 
reduced nonlinear activity at the small scales appears correct, the scaling of the observed 
spectrum is different than what was anticipated. 

A similar principle to reduce small-scale activity is used in the Leray model, ^^ which 
is based on Leray 's regular izat ion of the Navier-Stokes equations. ^° For the Leray model 



the momentum equations are the same as the Navier-Stokes equations with the exception 
that the advecting velocity is smoothed. Thus, they have the same form as equation (2) 
with the third term on the LHS set to zero and p" = p. However, in the Leray model it is 
the incompressibility of the unfiltered velocity which is enforced, diUi. It has been shown 
that this leads to significant problems when the model is used for wall-bounded fiows.^^ 

The reduced small-scale activity of the NS-a and Leray equations have led to the sugges- 
tion^^ to use these equations as turbulence models. This can be done either by working 
directly with the equations in which the unsmoothed velocity is the dependent variable, 
for example (2), or by rewriting these equations in terms of the smoothed velocity Mj. 
The latter is considered an LES methodology and gives rise to a governing equation that 
has the standard LES template. ^^ This is the approach that will be investigated here. 
Previous studies along these lines include the temporal transition of a mixing layer, ^^ 
the parameterization of mixing in a gyre,^^ and studies of decaying and forced box tur- 
bulence.^^ A conclusion that can be drawn from these studies is that the NS-a model 
captures the variability of the large, resolved scales better, ^^ or at least as well as,^^ the 
dynamic model. However, if the grid resolution is too coarse there will be a build-up of 
energy associated with the subgrid fiuctuations. This is particularly severe if the initial 
condition contains a broadband spectrum. ^^ These results are not surprising in light of the 
fact that while the NS-a model attenuates triad interactions associated with the forward 
energy transfer, it does not do so abruptly at a wavenumber corresponding to ka = 1/a. It 
is expected that a scale separation between the grid scale cut-off and a would be required 
to allow for this attenuation. If we consider that we still need to resolve the scales that 
eddies of scale 1/a transfer their energy to, this means kmax = 2/q;. Since the maximum 
wavenumber is also related to the grid spacing as kmax ~ tt/Zi this gives, a/h ~ 2/tt ~ 1. 
This is in agreement with the subgrid resolution suggested by Geurts and Holm,^^ which 
they determined through grid-refinement studies. 

While it is possible to define guidelines for how big a"^ should be in the case of homoge- 
neous, isotropic turbulence, it is not clear how to proceed in the more general situation. 
Since a^ is the only parameter in this model we expect its specification to be critical. 
Physically a^ can be interpreted either in the Lagrangian sense as a measure of rms par- 
ticle displacement, or in the Eulerian sense as a mixing length. Numerically, it can be 



interpreted as a filter width. ^^ The studies noted above all used a constant value for a^, 
which was taken to be a fraction of the domain size. It can be expected that there will 
be many situations where it may not be appropriate to maintain a constant value of a^, 
and in addition, where we may want a^ to refiect the anisotropy of the fiow. For example, 
near a solid wall we would expect the particle displacement and mixing length normal to 
the wall to decrease. In a similar manner, in LES we usually reduce the filter width and 
refine the grid in this region. 

As noted by Zhao and Mohseni*^^ there are two possible ways to proceed with the prob- 
lem of specifying a^ in an anisotropic fiow. One way would be to use the anisotropic 
NS-a equations^^'^^ to develop an equation for -Uj. A second approach would be to use 
the isotropic equations, but modify a^ to refiect the anisotropy of the fiow. Zhao and 
Mohseni followed the second approach, and formulated a dynamic procedure to specify 
(x^ . The model was tested a 'priori on channel fiow, where cP' was found to be constant 
away from the wall, decaying to zero such that a^ = at the wall. However, a posteriori 
results were less encouraging. "^° 

Outside of these studies by Zhao and Mohseni using the isotropic model, the author 
is not aware of any results in the literature where the NS-a model has been used in a 
numerical simulation of wall-bounded fiows of interest to engineers (although the Leray-a 
model has been studied, ^^ as have related symmetry-preserving methods^^). Given the 
importance of this problem, and the promise the model is showing for unbounded fiow, the 
objective of this study is to contribute to this picture by investigating using anisotropic 
the NS-a equations as a subgrid model. The anisotropic equations are comprised of a set 
of coupled PDEs governing momentum conservation and the particle displacement covari- 
ance.^^'^^ For the Eulerian-averaged equations investigated here, the governing equation 
for the displacement covariance is simply an advection equation, and it is not clear how 
such an equation would be initialized, or even if it should be treated in the prognostic 
sense. Instead of solving this equation, the initial approach taken here was to view o? as 
a smoothing scale which is based on the grid, and assess the performance of the model 
based on this definition. Unfortunately, difficulties were encountered with using a sim- 
ple mesh-based specification. For this reason an alternative definition of q\ was tested. 
It should be noted that results using the isotropic NS-a model are not presented here. 



We found this model generates excessive backscatter near the hd, leading to divergence. 
Similar problems have been found using the gradient or Clark model near walls, ^^'^^ and 
using the NS-a model at high values of (x^?'^ The outline of the paper is as follows. A 
description of the anisotropic subgrid model used is given in section 2. In this section 
condensed index notation is used to discuss the model for the sake of compactness. The 
numerical methods and treatment of the subgrid model are outlined in section 3. Results 
from using both a mesh-based definition of a^ and an alternative definition are presented 
in section 4. Concluding remarks are then given in section 5. 

2 Model Formulation 

The Eulerian-averaged equations from Holm^^ (section 12) are used as a starting point. 
In the development of these equations the instantaneous velocity is decomposed into a 
mean and a random fluctuation, and the averaging is applied at a fixed point. ^^ This is in 
contrast to the Lagrangian average, which is taken following a particle trajectory. We feel 
that the Eulerian average is more consistent with the manner in which the experimental 
and numerical data we are comparing with was collected. Differences between the Eulerian 
and Lagrangian averaged equations are discussed in Holm,^^ and alternative Lagrangian- 
averaged equations are given in Marsden.^^ The Eulerian-averaged equations are, 

diUi = 0, (5) 

dtUi + UjdjUi + UkdiUk = -diP + udkkUi - -di{^kCi)dkUmdiUm, (6) 

where P is a pressure-like variable, 

P = p- -UiUi - -{ikil)dkUmdlUm, (7) 

and with the following relationship between the smoothed and unsmoothed velocities, 

Ui={l-dk{{ikii)di))ui. (8) 

^ V ' 

H 

In (8) -u is a smoothed velocity, {S,k^i) is the smoothing scale, and the angle brackets, (■) 
denote an Eulerian average. For the isotropic model {^k^i) = ci^^ki- The last term on the 
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RHS of (6) arises in the derivation when the functional derivative is taken with respect to 
i^k^i) as is necessary to conserve momentum when {^k^i) is not constant. The momentum 
equation can also be written in momentum-conservation form as^ 

dtUi + UjdjUi = -diP + dj {{ikij)diUmdkUm) + vdkkUi. (9) 

While versions of the governing equations with the smoothed velocity as the dependent 
variable have appeared in the literature,^' ^^ the approach presented here is more familiar 
to the LES-community, and the purpose is to show that the subgrid stress rriij does not 
arise as an 'ad-hoc' modification to the isotropic model. To develop an equation with the 
smoothed velocity as the dependent variable we therefore follow the approach taken in 
Holm and Nadiga^^ and use the commutator between the substantial derivative and the 
smoothing operator. For example, we would like to have a substantial derivative written 
entirely in terms of the smoothed velocity. This is done by rewriting the advective terms 
in (9) as, 

dtUi + UjdjUi = [D/Dt, H]ui + H {dtUi + UjdjUi) . (10) 

Here [D/Dt, H] is the commutator between the material derivative and the Helmholtz 
operator, H from (8), [D/Dt,H]u = D/Dt{H{u)) ~ H{D/Dt{u)), where H{u) = u. Note 
that the substantial derivative is defined with the smoothed velocity, D/Dt = dt + Ujdj in 
keeping with the fact that the advecting velocity is smoothed. The momentum equation 
(9) can then be written as, 

dtUi + UjdjUi = H~^ {dip + dj {{^kQdiUmdkUm + i^OkkUi) - [D/Dt, H]ui) . (11) 

It was found the commutator can be expressed 



[dt + Ujdj, (1 - dk i{^k^i)di))] Ui = dj i{CkCi)dkUidiUj + {ijii)dkUidiUk) 

-9j {(dti^j^i) + Ukdki^j^i)) diUi) . 



(12) 



For constant, isotropic fluctuations, the first two terms on the RHS of (12) have the 
same form as the Leray model, ^^ although we used diUi = in the development of the 
commutator so this cannot be considered as the true Leray subgrid stress. For the NS-a 
model the last two terms on the RHS of (12) can in theory be neglected, because for the 
Eulerian-averaged equations each component of the particle displacement covariance is 



transported by the mean flow like a scalar, 

0. (13) 
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This arises in the derivation of the Eulerian-averaged equations when it is assumed that 
all of the fluctuations are contained in the Eulerian fleld.^^ 

In the context of LES modelling, this term represents the contribution to the subgrid 
stress due to the explicit change in fllter width. It is customary in the LES community 
to neglect these types of terms. For a recent discussion of this problem see van der Bos^^ 
(and references therein). As is the case with all subgrid terms, these terms will lead to 
energy transfer between resolved and subgrid modes, and neglecting these terms therefore 
changes the subgrid transfer dynamics. For the equation above, this is clearly seen by 
considering the isotropic case, where the last two terms on (12) can be written as, 

d ( ( Da^\ du 



dxk \\ Dt J dx 

The substantial derivative of a"^ can be seen to play the role of a variable eddy viscosity. 
It will dissipate energy when a^ increases along a flow path and backscatter energy when 
it decreases. This is exactly the method suggested in the literature^^ to model the com- 
mutation error in LES. The idea being that when a flow scale is advected into a region 
where the grid is coarser, it will go from being resolved to modelled, leading to dissipa- 
tion, and vice versa when the grid is reflned. In practice, whether or not they can be 
neglected depends on the magnitude of these terms relative to the other subgrid terms. 
In the one-dimensional case they can be neglected if 

1 da'^ 1 du 



a^ dx M (9x ' 

which means that the fllter fleld must be smoother than the flow fleld in the direction of 
advection. We found that the commutation term tends to be large when the advection is 
large, and does not have a signiflcant effect on the flow, as compared to the tilting term 
(which is UkdiUk in (2) and will become part of the dj{Cij) term in the following). 

The second simpliflcation will be to retain only the diagonal components of (^k^i)- This 



is equivalent to using the Helmholtz operator, 

H={1- dM^Qd^) - dy{{iyQdy) - 9, ( (^.C. ) 9, ) ) . (l4) 

This gives us a formulation similar to that derived using second order reconstruction 
methods, ^^ where the lack of off-diagonal terms arises when the three-dimensional filter 
is applied as the composition of three one-dimensional filters, L = h o h ° h^ where Ij 
{j = 1, 2, 3) represents a one dimensional filter in the Xj-direction. The benefit of such a 
simplification is that substantially reduces the cost and yields a subgrid model that is no 
more expensive than the isotropic version (which is itself about as costly as the dynamic 
modeP^). 

With the above assumptions and simplifications, the momentum equation can then be 
written 

dtUi + djUiUj = -dip + vdkkUi - H^^{djmij). (15) 



The subgrid stress is 



2 duiduj 2r dukdui 2 c du^dum 



where we have used al for {C,kC,k)- Following Geurts and Holm^^ the subgrid model can 
also be written as 

Here, Aij is the anisotropic gradient model (or Clark model), Aij + Bij is similar to a 
Leray model and the NS-a model is comprised of all three terms. 

2.1 Physical Interpretation of the Subgrid Term 

When written as a subgrid stress the effect of the rriij term is not easily interpreted 
physically. This term provides force to the momentum equation, and we will now rewrite 
it such that the form of this forcing function is clarified. Complementary discussions along 
these lines have been given in the literature,^' ^^ here we mention this explicitly because of 
its relevance to the following section. To keep things simple in this section we will assume 
isotropic fluctuations, and that a"^ is constant. We will also make use of the difference 



between the smoothed and unsmoothed velocity, 



wF = a^ 



dxl ' 



In the LES hterature this velocity would be called the subgrid fluctuation. In the NS-a 
literature it is referred to^^ as the 'Stokes velocity'. When a^ is constant the commutator 
can be written using the Stokes velocity as, 



[D/Dt, H] Ui = 2djAi 



u 



ST 



dui 



and the momentum equation can be written, 

dtUi + UjdjUi = —dip* + udkkUi — H~^ {2djAij — u^^ x uj) (19) 

where p* is a modified pressure. The term u^"^ x uj m (18) is called the vortex force. We 
can now see that the subgrid model is composed of two forcing terms. The first term, 
Aij^ is well known in the literature where it goes by many names, such as the Clark 
model, ^^ gradient model and Tensor-Diffusivity model. ^^'^^ It is a generic subgrid closure 
which can be derived by expanding the subgrid stress Tij in a Taylor series expansion and 
retaining terms up to O(A^), where A is the filter width. When the Helmholtz operator 
is approximated by a box filter, a^ ~ A^/24, where A is filter width, this term is identical 
to the model in the literature^^ 

24.. -2 2 9^i 9uj _ A^ duj duj 
*"' dxk dxk 12 dxk dxk 

This suggests an alternative (more approximate) way of deriving an equation for the 
smoothed velocity, which would be to start with the momentum equation, rewrite it with 
the Stokes vortex force on the RHS, apply a filter to the equations, and close the resulting 
Tij term with an explicitly filtered gradient model. 

The vortex force is what makes the NS-a model different from other approaches. To 
highlight how a vortex forcing term is fundamentally different than, for example, a 
Smagorinsky model, consider a simple two-dimensional mixing layer with u = tanh(|/). 
The Smagorinsky model will add a diffusion term to the momentum equation, with dif- 
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fusivity that is a function of the filtered strain-rate Sij, 

1/2 



UT = {CsAf (2S,,S,,') ~ {CsAf sech'y. 



The diffusivity will be highest at the middle of the mixing layer, and it is not surprising 
that such a model cannot be used for studies of mixing layer transition, where it damps 
out the small amplitude perturbations preventing transition. 

On the other hand, the vortex force would make its most significant contribution to 
the vertical (y) momentum equation, with a vertical forcing term 

C-rp 9 ^ ^ 

At the very early stages of transition this term would provide equal and opposite vertical 
forcing to the mixing layer, and therefore leaves the mixing layer unchanged. However, as 
soon as undulations in the layer appear the fiow is no longer symmetric and such a terms 
would serve to 'push' the mixing layer back and forth. Unlike the Smagorinsky model, 
the NS-a model was found to correctly capture the linear growth phase of a transitional 
mixing layer^^ 

3 Numerical Methods 

The governing equations for u are solved using the STREAM code.^^ This is a collocated 
finite-volume code which uses the SIMPLEC method to ensure mass conservation. The 
advection and diffusion terms are treated implicitly while the rriij term is treated using 
deferred correction. A second-order time stepping scheme is used with CFL ^ 1. The 
advection scheme used was the QUICK scheme. Note that this is very similar to methods 
used in earlier studies of this flow^^ and more recently in the evaluation of the dynamic 
mixed model on the same test case.^^ 

In the finite volume formulation the source term appears on the right-hand side of the 
momentum equation as 
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If we first calculate the rriij term at cell centers and then the divergence, this means the 
velocity gradients in the rriij terms would be computed at cell centers. It was found that 
doing this using second order central differences made the model even more susceptible to 
the numerical oscillations we will see in the next section. As an alternative, the method 
used here is to write the source term in a manner consistent with the other terms in the 
momentum equation 

H-^ [ ^ {mi,) dV - H-^ I m.jAj. (21) 

Jv (^^j Jcs 

The velocity gradients are now computed at control volume faces, which means the source 
term is computed from velocity differences between adjacent nodes. Thus the procedure 
is to first compute the subgrid force at all interior nodes, 

m.jAj, (22) 

cs 

and then find the filtered force Fi. For the present results this was done by solving the 
Helmholtz equation 

using a conjugate gradient solver. The solution of the Helmholtz equation requires bound- 
ary conditions for Fi, but in the discrete equation the boundary value of Fi will be mul- 
tiplied with a boundary value of al, which we know is zero from the boundary condition 
{Ck^i)nk = 0. Simulations were also carried out using a simple box filter, as has been done 
in previous studies. ^^ This was found be be much more efficient than Helmholtz inversion, 
with comparable results. 

4 Results 

4.1 Description of the Test Case 

The application of the NS-a model to a practical problem is studied here using a lid-driven 
cavity flow at a Reynolds number of 10,000, where the Reynolds number is based on the 
lid velocity and cavity length. The chosen cavity has a spanwise aspect ratio (SAR) of 
1, as shown in the schematic in Figure 1. The three-dimensional cavity flow contains a 
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variety of flow structures and is a challenging test case for a subgrid model due to the lack 
of homogeneous directions, the presence of both laminar and turbulent flow regions and 
the anisotropic nature of the flow. This cavity has been studied both experimentally^^ 
and numerically using LES^'^^ and DNS.^^ At this Reynolds number the distinguishing 
feature of the flow is the formation of two jets which separate off the downstream wall 
and impinge on the cavity bottom. While the experimental measurements reported a 
small inertial subrange near the cavity bottom, the DNS study reports that the flow does 
not actually become fully turbulent before it encounters the upstream wall. There are 
however, significant regions of both positive and negative turbulent energy production in 
the jet impingement regions, as is discussed in detail in the DNS paper. ^^ We would like 
to see to what extent the model investigated here can capture this. 

The mesh used is stretched in the x and y directions to capture the shear layers near 
the walls, but uniform in the spanwise since the relevant flow physics are not clustered 
near the endwalls, but distributed along the span. The parameters pertaining to the mesh 
sizes and stretching ratios are given in Table 1. The meshes are similar to those used in 
the study of Zang et al.^^ using the dynamic mixed model where similar numerical meth- 
ods were employed. 

To assess the time step we first compare our parameters to those used in the experiments. 
For the lid-driven cavity all quantities are non-dimensionalized by the cavity length (L) 
and lid velocity (U). The characteristic time scale is then L/U which can be written 
in terms of the Reynolds number as, -^ = -Re-^. Estimating the kinematic viscosity of 
water at room temperature as 1 x lO^^rri^/s and knowing the length of the cavity to be 
150mm gives L/U = 2.25s. A time step of 0.01 then corresponds to physical time of 
0.025s or physical frequency of 40Hz. The power spectra shown from the experiments for 
all cases have very little frequency content above IHz. Therefore, it was expected that 
the timestep chosen would be adequate. This was verified by the fact that simulation 
results showed very little difference when run at a timestep half as big. Rolling averages 
of the mean and rms velocities, as well as the total kinetic energy and dissipation were 
monitored. After a statistically steady state was achieved, statistics were collected over 
40L/f/. 
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4.2 Results with of, based on the mesh 

Since a^ is a smoothing scale we start with a simple definition based on the grid size 



tti 



C {hi) (24) 



where h^ is the grid spacing in the A;-direction and C is a constant denoting what fraction 
of the grid spacing to use. Because al can be related to the filter width, Aj^., of a box 
filter via al = A|/24/^ we choose C = 1/6, which corresponds to a filter width which 
is twice the grid size. It has been suggested that for adequate subgrid resolution in the 
isotropic model a"^ = K^ should be used for the NS-a model. ^^ In the present study the 
problems encountered were in laminar fiow regions, where the model should be inactive, 
and the question of subgrid resolution was not addressed in detail. However, in some 
cases simulations were run on different meshes to verify the sensitivity of the results to 
the observed trends. In all cases the value of C was adjusted so that the physical value 
of a\ was approximately the same as on the coarse mesh. Simulations were also done 
with the isotropic model with a^ proportional to the volume, but as mentioned in the 
introduction, this model was found to be generate excessive backscatter near the cavity 
lid, leading to divergence. 

It was found that there was a very persistent problem when a^ was based on the mesh. 
This was that the downstream wall jet was pushed too far out from the wall, as shown 
in Figure 2. This was observed on both coarse and refined meshes, over a range of a\ 
values, and also when a box filter was used instead of a Helmholtz operator. It was also 
seen when the isotropic version of the model (with a^ based on the grid volume) was used. 

Because the wall jet is pushed too far out from the downstream wall (in the x-direction) 
we now look at the contribution of the subgrid force to the -u-momentum equation. Using 
the vortex force. 

In the wall jet region we found the vortex forcing term is much larger than the Aij term. 
In this region, oj^ ^ <^y and d^ ^ dy, dz and it was expected that the problem was coming 
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from following component of the vortex force, 

Given that the vorticity field is unsteady this will be an unsteady forcing term which 
could cause the wall jet to oscillate back and forth, leading to high fluctuation levels. 
Depending on the balance between the positive and negative forcing, it is also possible 
that this could lead to the jet being pushed too far out from the wall in the mean. For 
the anisotropic model this hypothesis could be easily tested by turning off a^. To our 
surprise this did not help the situation. Instead, it was turning off a1 which solved this 
problem. 

A closer examination of the flow fields corresponding to the a1 ^ and a^ = cases 
showed that the main difference between the two is the appearance of streamwise (here 
vertical) vortices (ujy) in the downstream walljet region when a^ ^ 0, as shown in Figure 
3. These vortices do not appear when a model is not used and appear to be a numerical 
artifact of the NS-a model. These presence of these vortices can be understood if there is 
significant modulation of the velocity in the spanwise direction, as for example could be 
caused by spanwise numerical oscillations. Recall that the vortex force in the momentum 
equations appears as advection and stretching/tilting terms in the vorticity equations. In 
particular the stretching/tilting term in the vertical vorticity equation is, 

QyST ^ Q^ST ^ Q^ST 
^x^ \- ^y^7\ 1" ^z^T, • 

OX oy oz 

Oscillations in the spanwise direction would lead to the generation of streamwise vorticity, 
tilted into the vertical by the first term. It was found that the fiat velocity profile seen in 
Figure 2 developed slowly over time, indicating that the long-time average effect of the 
vertical vortices is a mixing out of the velocity profile. Similar problems were observed in 
cases where the numerical oscillations were not as visually obvious (at a lower Reynolds 
number of Re = 3,200), again showing that since the effect builds up over time, small 
oscillations can have a significant impact. 

There is an interesting analogy between the behavior seen here, and the true physical 
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behavior of the NS-a equations. It has been pointed out^^ that these equations are re- 
lated to the Craik-Leibovich equations/ which are used to account for the long-time 
averaged effect of surface waves on the background current. The effect of surface waves 
is to create a relative velocity between a fluid particle (Lagrangian) and the background 
current (Eulerian). This relative velocity, which is called the Stokes drift velocity, then 
tilts vertical vorticity into the streamwise direction to create streamwise vortices (Lang- 
muir cells) which transport momentum perpendicular to the free surface and flatten the 
velocity profile below the surface, leading to a mixed layer, ^^ a schematic is shown in 
Figure 4. Thus, this result is not erroneous in the sense that it is a real solution to the 
given equations in the presence of small scale spanwise oscillations. However, the problem 
is that the oscillations are not coming from something physical, such as surface waves, 
but from an unwanted numerical effect. 

To verify that this is not just a boundary problem a simple test was done with a laminar 
Couette flow and an a^ discontinuity was introduced in the middle. Once again, the 
model generated a force due to the C22 term, which was balanced by a pressure gradi- 
ent. This was tested with varying subgrid resolutions (i.e. keeping the physical size of ai 
constant but refining the mesh) and was found to be relatively insensitive to the resolution. 

Although this is not a boundary problem per se, it can be readily verified that if both 
^ ^ and dyol -^ at the lid, we would not have this problem. From the relationship 



a 



between the smoothed and unsmoothed velocity fields, 

dy dy ^ dy'^ ' 

we can see this corresponds to both fields satisfying the same boundary condition, which 
here would be the no-slip condition. Unfortunately when a^ is based on the mesh and 
the mesh is uniform it is impossible to satisfy both a^ — » and dya^ — * 0. Note that 
the same problem arises for the isotropic version of the model, which was also verified 
numerically. 
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4.3 An alternative definition of al 

As an alternative to having al based strictly on the mesh spacing, we have used a defi- 
nition of al which incorporates the properties of the resolved flow. There are a number 
of ways one could do this, and the definition given here is one which we found to work well. 

To gain further insight in how to define a^. we can go back to Taylor's paper^^ on turbulent 
diffusion since, in theory, a^ is a measure of mean-squared particle displacement. The 
relationship between the NS-a model and Taylor's work on turbulent diffusion has been 
discussed earlier.^''' First of all, Taylor showed that if the averaging time for the particle 
motion (T) is long relative to the time over which the particle takes a step (r) the scaling 
of the mean-squared particle displacement will be, [X^] ~ [f]^T^, where [v]"^ is a measure 
of the particle velocity and [■] denotes an ensemble average. In contrast if the averaging 
time is short relative to the step time T r^ r the scaling will be, [X^] ~ [v^'Tt. In the 
development of the NS-a model it is assumed that there is a separation of scales. ^^ Thus 
we will apply the former scaling here and use a time scale based on the velocity gradient 
tensor T^ ~ QijQij where gij = djUi. To form a velocity scale we again follow Taylor^^ 
who pointed out that in considering the dispersion of a particle due to turbulent motion 
it is not the kinetic energy of the particle f ^ that is relevant, but the number of times it 
changes direction. In one dimension this can be captured by [d^v) or {dfv) . In the more 
general case a second-order structure function could be used. In the anisotropic case this 
would be,^^ 



-"5 



3 X , x2/3 



1 / A \' 

F2(x,A,t) = -^[||u(x,t)-u(x + Aa;,e„t)||2-||u(x,t)-u(x-Aa;,e„t)||2] (^j 

(27) 
Here ej denotes a unit vector and A is a length scale based on the grid volume as 
A = (/ii/i2^3) • For homogeneous, isotropic turbulence this is similar to using the 
turbulent kinetic energy to estimate [v]'^ since in that case there is a simple relationship 
between the second order structure function and the energy spectral density (see Batch- 
elorp. 120^). 

Putting the velocity and time scales together we would then arrive at the following defi- 
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nition for a^, 



2 i'2(x,A,t) 



ai = ^^V^. (28) 

In practice F2 is computed using the six closest neighbors to a given mesh point. ^ This 
means such a definition of a^ would reduce to the wall normal spacing in a wall-bounded 
flow, which will result in little improvement over the simple grid-based definition. This 
problem can be anticipated because in a wall bounded fiow, for example a channel fiow 
with du/dy as the shear, the velocity fiuctuation associated with uiy + Ay) — u{y) is not 
fully turbulent, and should not be included in the computation of F2. This problem has 
been discussed in the literature in applications of the structure function model to chan- 
nel and boundary layer fiows.^ In this case the problem was resolved by not including 
uiy + Ay) —u{y) in the calculation of F2. In the more complex situation other strategies, 
such as high pass filtering, are often used.^ 

The definition which was found to work well instead was, 

al = max [{5^u)\ {5yu)\ {5,u)^] T^ (29) 

[{5.v)\ {6yv)\ {6,vf] T' (30) 



ay = max "^ --^ ' ^ -^ /^ „.\ 1 t^ 



al = max [{6,w)\ ((5»^ {6,w)^] T^ (31) 



where again T^ is {gijgij)~^ and the S symbol denotes a velocity increment. In practice 
this can computed as the velocity difference between adjacent mesh points. Whereas a 
structure function is based on the velocity difference in a given direction and tells us about 
energy contained in eddies of a given size, this definition tells us about the energy in the 
horizontal, vertical and spanwise velocity fiuctuations. The question then arises as to 
which is more appropriate. The definition given above was based on heuristic reasoning. 
If a blob of fiuid is experiencing an oscillating shear force, it would be the du/dy shear 
which would cause it to move back and forth in the horizontal direction, while the dv/dx 
shear would cause it move back and forth in the vertical direction. Thus it was reasoned 
that al should be related to 5uy and not 5vx- 



4.4 Results from the alternative definition 

We now look at the performance of the model with the alternative definition of al given in 
equations (29)-(31). For comparison, results are also shown for the case where no subgrid 
model is used. There are several ways the performance of a subgrid model can be assessed. 
We start by looking at how well the mean flow is captured, which is reflected in the wall jet 
structure. Recall that the flow should split into two wall jets, which impinge on the cavity 
bottom. We can see in Figure 5 that when a model is not used the flow does not split into 
two jets, and that this situation is corrected when we used the NS-a model. We found 
that even on the coarse mesh of (32)^ the NS-a model with the alternative definition of a^ 
can correctly produce the splitting into two wall jets. However, the energy spectra at such 
a coarse resolution did not exhibit a k~^^ slope, so no results from this test case are shown. 

The mean flow, rms and shear stress profiles are shown in Figure 6 for the 48^ mesh 
and in Figure 7 for the 64^ mesh. In Figure 6 we also show the profile from using the 
mesh-based a]:. It can be clearly seen the flow-dependent definition is necessary to obtain 
the correct mean flow profile. It can also be clearly seen that the new model does a good 
job of capturing the velocity fluctuations near the lid and in the downstream wall jet 
region, and that the shear stress proflles are in excellent agreement with the experimental 
data.^^ In contrast, without the model (solid line) the fluctuations are too low, and the 
shear stress is underpredicted. For the flner mesh results shown in Figure 7 the differences 
with and without the model are small, indicating that as a^ — > the simulation moves 
towards a DNS as it should. 

The highly inhomogeneous and anisotropic nature of lid driven cavity flow has been was 
well documented in the DNS and LES studies of Leriche and Gavrilakis^^ and Bouffanais 
and Deville.^ One measure of anisotropy they used is the ratios of the volume- averaged 
contributions of the mean velocity components to the kinetic energy. In the present study 
it was found the ratio jy{u)'^ dV : jy{v)'^ dV : jy{w)'^ dV was 1 : 1.23 : 118 without the 
model as compared to 1 : 1.21 : 60 with the model, both on the 64'^ mesh. The model 
compares much more favorably with the DNS study which reported 1 : 1.22 : 50. This can 
be expected from the stronger impingement of the wall jet when the model is used, and 
the resulting momentum transfer into the spanwise direction. The stronger impingement 
is very evident if we look at the contours of the production term, P22 = —{v'v')dy{v). The 
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contours shown in Figure 8 are in good qualitative agreement with the DNS study^^ (Figure 
14). 

Since the flow in the downstream wall jet region is characterized by positive and neg- 
ative turbulent energy production^^ we expect the contribution of the subgrid model to 
the resolved flow energy equation to exhibit positive and negative values in this region 
also. The contribution of the subgrid stress to the resolved flow energy equation is 

^ drriij d dui 

Ui^ — = 7^— [Uimij) - rriij-—. 

OXj OXj OXj 

The first term on the RHS is the transport due to the resolved flow while the second 
is a source/sink term, usually referred to as the SGS dissipation term. Since it can be 
both positive or negative, we prefer to call it the SGS transfer term, as it is responsible 
for the energy transfer between the resolved and subgrid modes (there is an equal and 
opposite term in the subgrid-scale energy equation'^'^). In our method we do not compute 
rriij explicitly, but rather the volume-integrated subgrid force, 

Jv (y^j 

This means we cannot split the energy transfer into these two contributions but instead 
plot the total SGS contribution, UiFi divided by the control volume. Contour plots of 
this term on a plane near the cavity bottom are shown in Figure 9. It can be seen there 
are both negative and positive contributions, and that the impingement points are asso- 
ciated with the energy transfer from the resolved flow, while the spreading is associated 
with energy transfer to the resolved flow. This is in good agreement with the DNS which 
found both positive and negative turbulent kinetic energy production terms in this region. 

To compare the current definition of a^ given in (29)- (31) with the mesh-based defi- 
nition from equation (24), plots of a\/h\ are shown in Figure 10. We can see that 
c^l/hy is high in the jet impingement region, while a^/hl and al/h^ reflect the spreading 
of the jet on cavity bottom, and the impingement on the upstream wall. Considering 
that the relationship between the unsmoothed and smoothed velocity in Fourier space is 
Ui{k) = (1 + a^k'^)ui{k) and the maximum resolvable wavenumber is A; ~ n/h we can also 
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look at this as the range of (ak) values. When (ak) = the model is inactive, while in 
the turbulent regions we expect (ak) ~ 1. This is reflected in the plots shown in Figure 
10. 



The actual force experienced by the flow due to the subgrid model is also of interest. 
In Figure 11 we plot the subgrid force contribution to the a;— momentum equation, which 
can be compared to the mesh based deflnition discussed earlier. It can be seen that the 
high source terms near the lid and in the downstream wall jet region are eliminated when 
the flow dependent version of a]: is used, and instead the flow is active in the turbulent 
regions near the cavity bottom. 



5 Conclusions 

An anisotropic version of the NS-a subgrid model (where u is the dependent variable) was 
developed starting from the anisotropic Eulerian-averaged equations given by Holm^^ in a 
manner that should be familiar to the LES community. While simpliflcations were made 
(Section 2), this work here still represents (to the author's knowledge) the flrst application 
of the anisotropic NS-a equations as a subgrid model in the context of LES, including the 
solution of a wall-bounded flow. Because the isotropic equations are showing promise in 
unbounded flows, we view this as a flrst step towards a more general application of the 
model to complex flows. The full anisotropic subgrid model should be tested against the 
one used here (with only the diagonal al) using a more appropriate test case in a separate 
study. 

The model was found to be sensitive to abrupt changes in a^. This is not surprising 
since a"^ is supposed to be a smoothing parameter, and abrupt changes are hardly physi- 
cal. However, if a^ = on the solid boundary is to be enforced, it was found this can be a 
problem. For the three-dimensional cavity flow this problem manifested itself in the form 
of oscillations in the spanwise velocity fleld and in the appearance of small-scale vertical 
vorticity. This vorticity can be understood as being due to the tilting of the spanwise vor- 
ticity from the Stokes-vortex force, an effect here which is numerical rather than physical. 
To overcome this problem an alternative deflnition of al was proposed which is not based 
solely on the mesh spacing. This deflnition worked very well in capturing the wall-jet 
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splitting, rms and shear stress profiles, and was also found to predict both forward energy 
transfer and backscatter in the jet impingement regions in qualitative agreement with the 
discussion given in Leriche and Gavilakis.^^ The alternative definition allows us to use the 
model in a complex flow situation that presents a significant challenge to most subgrid 
models. For the lid-driven cavity it was important that the model remain inactive in the 
laminar flow regions, which was not possible when a^ was based on the mesh. While this 
was not a problem in the mixing layer study carried out by Geurts and Holm,^^ it should 
be noted their problem did not have solid boundaries, and was relatively symmetric in 
the early stages of transition. 

Lastly it should be mentioned that simulations were also done with the C^ term turned 
off, which is similar to using a Leray model. It was found in these cases that there was 
no benefit to using the model, and in some cases the model tended to damp the small 
scale activity strongly. This is in agreement with recent results^^ which indicate the Leray 
model reduces the effective Reynolds number of the flow. The tilting term, UkdiUk, which 
combines with the modifled pressure gradient to form the dj{Cij) term in the model, is 
the unique feature of the NS-a model. The role of this term is presently being investi- 
gated in turbulent channel flows. It is hoped the channel flow cases will also delineate the 
near-wall behavior of the model further. 
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Figure 1: Sketch of the hd-driven cavity flow. 
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Table 1: Mesh parameters 



2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 




Figure 2: Mean flow profiles on the midplane for the (64)^ mesh showing the wall jet is 
pushed out too far from the downstream wall when mesh-based al is used. Solid line, no 
model; dashed line, NS-a model with al based on the mesh. Symbols are experimental 
data.^^ 
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(a) a1 based on the mesh 




(b) al - 

Figure 3: Vertical vorticity ujy near the downstream wall demonstrating the small-scale 
vorticity found due to the al discontinuity. The plane is at a height oi y = 0.6. 
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Figure 4: A spanwise perturbation generates vertical vorticity. Opposing torques t are 
created from the u^'^ x uj term in the momentum equation. These create a convergence 
zone and vertical vorticity. In the vorticity equation, horizontal vorticity (uj^) is tilted by 
v^'^ into the vertical direction. After the sketch in Leibovich/^ rotated 90°. 
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Figure 5: Contours of (v) looking down at the cavity bottom showing that the wall jet 
correctly splits into two when the model is used, but does not split without the model. 
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(a) Mean flow 




rms 
(b) rms profiles 




(c) {u'v') profiles 

Figure 6: Mean flow, rms and {u'v') profiles on the midplane for the 48^ mesh. Sohd fine 
is no model, dotted line (first plot only) is NS-a with the mesh-based definition of al and 
dashed line is NS-a with alternative definition of al- Symbols are experimental data.^^ 
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Figure 7: Mean flow, rms and {u'v') proflles on the midplane for the 64^ mesh. Sohd hne 
is no model, dashed hne is with alternative deflnition of a^. Symbols are experimental 
data.^^ 




Figure 8: 
0.045. 



P22 contours on the y = 0.03 plane for the 64^ mesh, levels between —0.015 and 
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Figure 9: Energy transfer term UiF^ on the y = 0.02 plane for the 64^ mesh. 
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(a) {ay/hy) on the z = 0.3 plane 




(b) {al/hl) on the y = 0.01 plane 




(c) (al/hl) on the y = 0.01 plane 

Figure 10: Contour plots of ot^Jh'^j highlighting the wall jet impingement and spreading 
regions for the 64'^ mesh. 
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(b) alternative definition of o?^ 

Figure 11: Subgrid force to the x— momentum equation on the z = 0.3 plane for the 64^ 
mesh. With al based on the grid the force is high in the laminar regions (near the lid 
and downstream wall), whereas with the alternative definition (equations (29)-(31)) the 
force is high only in the turbulent regions. 



29 



References 

[1] G.K. Batchelor. The Theory of Homogeneous Turbulence. Cambridge University 
Press, 1972. 

[2] R. Bouffanais and M. O. Deville. Large-eddy simulation of the flow in a lid-driven 
cavity. Physics of Fluids, 19:055108, 2007. 

[3] S. Chen, C. Foias, D.D. Holm, E Olson, E.S. Titi, and S. Wynne. Camassa-Holm 
equations as a closure model for turbulent channel and pipe flow. Physical Review 
Letters, 81:5338-5341, 1998. 

[4] S. Chen, C. Foias, D.D. Holm, E. Olson, E.S. Titi, and S. Wynne. The Camassa-Holm 
equations and turbulence. Physica D, 133:49-65, 1999. 

[5] S. Chen, C. Foias, D.D. Holm, E. Olson, E.S. Titi, and S. Wynne. A connection 
between the Camassa-Holm equations and turbulent flows in pipes and channels. 
Physics of Fluids, 11:2343-2353, 1999. 

[6] S. Chen, D.D. Holm, L.G. Margolin, and R. Zhang. Direct numerical simulations of 
the Navier-Stokes alpha model. Physica D, 133:66-83, 1999. 

[7] A. D.D. Craik and S. Leibovich. A rational model for Langmuir circulation. Journal 
of Fluid Mechanics, 73:401-426, 1976. 

[8] J. A. Domaradzki and D.D. Holm. Navier-Stokes alpha model: LES equations with 
nonlinear dispersion. In B.J. Geurts, editor. Modern Simulation Strategies for Tur- 
bulent Flow, chapter 6. R.T. Edwards, Inc., 2001. 

[9] F. Ducros P. Comte and M. Lesieur. Large-eddy simulation of transition to turbulence 
in a boundary layer developing spatially over a flat plate. Journal of Fluid Mechanics, 
326:1-36, 1996. 

[10] C. Foias, D.D. Holm, and E.S. Titi. The Navier-Stokes- alpha model of fluid turbu- 
lence. Physica D, 152-153:505-519, 2001. 

[11] C.J. Freitas and R.L. Street. Non-linear transient phenomena in a complex recircu- 
lating flow: A numerical investigation. International Journal for Numerical Methods 
m Fluids, 8:769-802, 1988. 

30 



[12] B.J. Geurts. Elements of Direct and Large-Eddy Simulation. R.T.Edwards, 2003. 

[13] B.J. Geurts and D.D. Holm. Leray and LANS-alpha modelling of turbulent mixing. 

Journal of Turbulence, 7(10): 1-33, 2006. 

[14] J. Graham, D. Holm, P. Mininni, and A. Pouquet. Three regularization models of 
the Navier-Stokes equations. Physics of Fluids, 20:035107, 2008. 

[15] K. Hanjalic. Will RANS survive LES: a view of perspectives. Journal of Fluids 
Engmeermg, 127:831-839, 2005. 

[16] D.D. Holm. Fluctuation effects on 3D Lagrangian mean and Eulerian mean fluid 
motion. Physica D, 133:215-269, 1999. 

[17] D.D. Holm. Taylor's hypothesis, Hamilton's principle and the LANS-alpha model 
for computing turbulence. Los Alamos Science, (29):172-180, 2005. 

[18] D.D. Holm and B.T. Nadiga. Modeling mesoscale turbulence in the barotropic 
double-gyre circulation. Journal of Physical Oceanography, 33:2355-2366, 2003. 

[19] S. Leibovich. The form and dynamics of Langmuir circulations. Annual Review of 
Fluid Mechanics, 15:715-724, 1983. 

[20] J. Leray. Sur les movements d'un fiuide visqueux remplissant I'espace. Acta Mathe- 
matica, 63:193-248, 1934. 

[21] E. Leriche and S. Gavrilakis. Direct numerical simulations of the flow in a lid-driven 
cubical cavity. Physics of Fluids, 12:1363, 2000. 

[22] M. Lesieur and O. Metais. New trends in large-eddy simulations of turbulence. 
Annual Review of Fluid Mechanics, 28:45-82, 1996. 

[23] F.S. Lien and M.A. Leschziner. A general non-orthogonal collocated FV algorithm for 
turbulent flow at all speeds incorporating second moment closure. Part 1: Computa- 
tional implementation. Computer Methods for Applied Mechanics and Engineering, 
114:123-148, 1994. 

[24] J.E. Marsden and S. Shkoller. The Anisotropic Lagrangian Averaged Euler and 
Navier-Stokes equations. Archives of Rational Mech. Analysis, 66:27-46, 2003. 

31 



[25] J.C. McWilliams, P.P. Sullivan, and CH. Moeng. Langmuir turbulence in the ocean. 
Journal of Fluid Mechanics, 334:1-30, 1997. 

[26] K. Mohseni, B. Kosovic, S. Shkoller, and J.E. Marsden. Numerical simulations of the 
Lagrangian Averaged Navier-Stokes equations for homogeneous isotropic turbulence. 
Physics of Fluids, 15(2):524-544, 2003. 

[27] M.R. Petersen, M.W. Hecht, and B.A. Wingate. Efficient form of the LANS-alpha 
turbulence model in a primitive equation ocean model. Journal of Computational 
Physics (accepted), 2008. 

[28] U. Piomelli, W.H. Cabot, P. Moin, and S. Lee. Subgrid-scale backscatter in turbulent 
and transitional flows. Physics of Fluids A, 7(3):1766-1771, 2001. 

[29] A.K. Prasad and J.R. Koseff. Reynolds number and end-wall effects on a lid-driven 
cavity flow. Physics of Fluids A, 1(2):208-218, 1988. 

[30] P. Sagaut. Large Eddy Simulation for Incompressible Flows. Springer- Verlag, 2002. 

[31] G.l. Taylor. Diffusion by continuous movements. Proceedings of the London Mathe- 
matical Society, 20:196-212, 1922. 

[32] F. van der Bos and B.J. Geurts. Commutator errors in the filtering approach to 
large-eddy simulation. Physics of Fluids, 17:035108, 2005. 

[33] M. van Reeuwijk. Direction simulation and regularization modeling of turbulent ther- 
mal convection. PhD thesis. Delft University of Technology, 2007. 

[34] M. van Reeuwijk, H.J.J. Jonker, and K. Hanjalic. Incompressibility of the Leray- 
alpha model for wall-bounded flows. Physics of Fluids, 18:018103, 2006. 

[35] Verstappen2007. Symmetry-preserving regularization modeling of turbulent channel 
flow. In Turbulent Boundary Layers, June 2007. 

[36] B. Vreman. Comment on Tnapplicability of the dynamic Clark model to the 
large eddy simulation of incompressible turbulent channel flows. Physics of Fluids, 
16(2):L29, 2004. 



32 



[37] B. Vreman, B.J. Geurts, and H. Kuerten. Large-eddy simulation of the temporal 
mixing layer using the Clark model. Theoretical and Computational Fluid Dynamics, 
8:309-324, 1996. 

[38] G.S. Winckelmans, O. Wray, A. A. Vasilyev, and H. Jeanmart. Explicit-filtering 
large-eddy simulation using the tensor-diffusivity model supplemented by a dynamic 
Smagorinsky term. Physics of Fluids, 13(5):1385-1403, 2001. 

[39] Y. Zang, R.L. Street, and J.R. Koseff. A dynamic mixed subgrid-scale model and its 
application to turbulent recirculating flows. Physics of Fluids A, 5(12):3186-3196, 
1993. 

[40] H. Zhao and K. Mohseni. Anisotropic turbulent flow simulations using the 
Lagrangian- Averaged Navier-Stokes alpha equation. In Proceedings of the 15th AIAA 
Fluid Dynamics conference and Exhibit, June 2005. 

[41] H. Zhao and K. Mohseni. A dynamic model for the Lagrangian Averaged Navier- 
Stokes a equations. Physics of Fluids, 17:075106, 2005. 



33 



